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Abstract. Recently, there has been plenty of work in designing and fabricating materials with 
an effective negative refractive index. Veselago realized that a slab of material with a refractive 
index of —1 would act as a lens. Pendry suggested that the Veselago lens would act as a superlens, 
providing a perfect image of an object in contrast to conventional lenses which are only able to focus 
a point source to an image having a diameter of the order of the wavelength of the incident field. 

Recent work has shown that similar focusing effects can be obtained with certain slabs of "con- 
ventional" periodic composite materials: photonic crystals. The present work seeks to answer the 
question of what periodic dielectric composite medium (described by dielectric coefficient with pos- 
itive real part) gives an optimal image of a point source. An optimization problem is formulated 
and it is shown that a solution exists provided the medium has small absorption. Solutions are 
characterized by an adjoint-state gradient condition, and several numerical examples illustrate both 
the plausibility of this design approach, and the possibility of obtaining smaller image spot sizes than 
with typical photonic crystals. 
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1. Introduction. Recently, there has been a renewed and avid interest in study- 
ing a class of materials known as the left-handed materials (LHMs) . These materials 
have simultaneously negative real parts of dielectric permittivity p and magnetic per- 
meability /i, so that their refractive index is negative. The properties of such meterials 
were investigated first by Veselago in 1967 [20] . As shown by Veselago, LHMs exhibit 
some peculiar electromagnetic properties such as negative index of refraction and wave 
vector, k, and Poynting vector, S, having opposite directions. Veselago realized that 
a slab of LHM would act as a lens. 

Due to the absence of naturally ocurring materials possessing both negative per- 
mittivity and negative permeability, Veselago's predictions did not receive much at- 
tention until recently, when a material with both negative permittivity and negative 
permeability at microwave frequencies was built [12] • Subsequently, the properties of 
LHMs were analyzed by many authors. 

According to Abbe's diffraction limit, conventional lenses based on positive index 
materials with curved surfaces are not able to resolve an object's fine details that are 
smaller than half of the light wavelength A. The limitation occurs because the waves 
with transverse wave numbers larger than 27m/ A, which carry information about the 
fine sub-A details of the object, decay exponentially in free space. In a negative index 
material slab, however, the evanescent wave components can grow exponentially and 
thus compensate for the exponential decay. Therefore, under ideal conditions, all 
Fourier components from the object can be recovered at the image plane producing a 
resolution far below the diffraction limit [17 . 

Milton et al. proved superlensing in the quasistatic regime (where the wavelength 
is much larger than the object), and discused limitations of superlenses in this regime 
due to anomalous localized resonance. If the source being imaged responds to an 
applied field, it must lie outside the resonant regions to be successfully imaged [14]. 
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In the electrostatic limit, the magnetic and electric fields decouple, and the re- 
quirement for super lensing of transverse-magnetic waves is reduced to only p — —ph, 
where ph is the permittivity of the host medium interfacing the lens [17 . An example 
of such near field superlens is a slab of silver in air illuminated at its surface plas- 
mon resonance (where p ~ — !)■ Experiments with silver slabs have already shown 
rapid growth of evanescent waves |10| . submicron imaging 13 , and imaging beyond 
the diffraction limit [6 . A major draw-back of such near-field superlenses based on 
bulk metals is that they can operate only at a single frequency u) satisfying the lens 
condition p{lo) = —ph- Shalaev et al. proposed a "tunable" near-field superlens made 
of metal-diclcctric composites that can operate at any desired visible or near-infrared 
wavelength with the frequency controlled by the metal filling factor of the composite 
(here the inhomogeneities are assumed to be much smaller than the wavelength) [T]. 

It was shown by Efros et al. that a two dimensional photonic crystal made 
from a non-magnetic dielectric has negative values of both the electric permittivity 
and the magnetic permeability in some frequency range 4 , and the photonic crystal 
behaves like a real LHM with respect to the propagating modes only. If amplification 
of evanescents waves occurs, it is due to some other reason; for example, excitation 
of surface waves by the evanescent waves. Such an amplification may provide an 
improvement of the image in the near-field region, but it does not affect the image 
near the far-field focal point [5]. 

The physical principles that allow negative refraction in photonic crystals arise 
from the dispersion characteristics of wave propagation in periodic media and are 
very different from those in LHMs. They also do not require both negative electric 
permittivity and magnetic permeability [3 [16] . The negative refraction of beams can 
be described by analyzing the equifrequency surface of the band structures O [HI [H] . 
If the constant-frequency contour is everywhere convex, an incoming plane wave from 
air will couple to a single mode that propagates into the crystal on the negative side 
of the boundary, and thus negative refraction in the first band is realized. Luo et 
al. have shown all-angle negative refraction could be achieved at the lowest band 
of two-dimensional photonic crystals in the case of S • k > [11]. Such all-angle 
refraction is essential for superlens application. The photonic crystal not only focuses 
all propagating waves without limitation of finite aperture, but also amplifies at least 
some evanescent waves, and the unconventional imaging effects are due to the presence 
of additional near-field light. A perfect lens, made of left-handed materials, focuses all 
propagating waves and all evanescent waves. The important difference for superlensing 
with a photonic crystals is that only finite number of evanescent waves is amplified. 
This is a consequence of Bragg scattering of light to leaky photon modes [H]. The 
resolution of a photonic-crystal superlens at a single frequency is only limited by its 
surface period instead of the wavelength [12]. 

More recently Huang et al. proposed an alternative approach to all-angle negative 
refraction in two-dimensional photonic crystals. By applying appropriate modifica- 
tions with surface grating to the flat photonic lens, he is able to focus large and/or 
far way objects 0. 

Inspired by the current research in structures that produce sub- wavelength focus- 
ing, we use derivative-based minimization techniques to produce structures that will 
provide sub-wavelength focus with non-magnetic materials and without the need for 
negative permittivities. Rather than restricting to designs based on photonic crys- 
tal structures, we allow as admissible any periodic composite structure (with fixed 
period) whose refractive index is bounded above and below by fixed constants. And 
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rather than performing parametric optimization over a smaU number of variables de- 
scribing the structure, we use techniques of "topology optimization" in which material 
distribution is completely arbitrary. We are able to obtain structured that focus a 
point source that is far away from the lens, and also we can obtain structures that 
give an image at a chosen distance from the lens. Since structures incorporating grat- 
ings are included in our admissible class, such designs will naturally arise through the 
optimization process if they produce the best possible image. 

For simplicity, only the case of "two-dimensional" structures in i?-parallel polar- 
ization is considered. The ideas here should extend to the other polarization case and 
the full three-dimensional problem, although there are some technical hurdles. 

The paper proceeds as follows. In Section 2 we describe the model problem and 
review a variational formulation of the Helmholtz equation in a periodic geometry. 
The inclusion of a small amount of energy absorption in the medium allows a uniform 
upper bound on the norm of the electric field, independent of the particular admis- 
sible structure (and thereby preventing resonances). In Section 3 we present the 
optimization problem and derive the optimality conditions. In Section 4 we depict 
the numerical experiments and show structures that have produced sub-wavelength 
focus. 

2. Model problem. In this paper we consider time-harmonic electromagnetic 
wave propagation through nonmagnetic {fi — 1) heterogeneous media for which the 
dielectric coefficient is constant in one direction, i.e. e(x,y,z) — p{x,y). Assuming 
that the electric field vector E ~ (0, 0, u), Maxwell's equations reduce to the Helmholtz 
equation 



where oj represents the frequency, and p G L°^{M.'^) is the dielectric coefficient. 

2.1. Periodic structure. Assume that the dielectric coefficient p{x,y) is peri- 
odic in the x variable 



Taking the period to be 2tt imposes no loss of generality since any other period can 
be obtained by rescaling u. 

Assume that the regions {y > 0}, and {y < —b} are homogeneous, for some fixed 
constant 5 > 0. In particular, assume for y > and y < —6, that p{x,y) — 1. The 
slab —b < y < may contain inhomogeneous material. 

Suppose a point source is placed above the slab at the point (0, h), which generates 
the incident field Ui{x, y) — H'^\ijjr), where r = ^ x"^ -\- {h — y)^ is the distance from 
the source, and is the Hankel function. For y < h, we have the representation 



where (3{£_) — uj^ — whenever the argument is positive, and — i^/S,"^ — ui'^ 
otherwise p^. It follows that 



Au + Lu'^pu = 0, inR2, 



(2.1) 



p(x, y) = p{x + 27r, y), for all (x, y) € M^. 
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and 



We can rewrite 



inx 



and 



/•1/2 

g(x)^ gJx)e'"^da, where 5^(2;) = — V e*''("+")''e 

^-1/2 TT f:^ 



'i/3(n+a)/i ina; 



{fa may fail to converge, but only at the isolated values of a for which (n + a)^ = 
for some n. always converges due to the exponential decay in n.) 

Above the slab y > 0, we separate the solution u to (|2.ip into the incident and 
scattered field: u — Ui + Ug- The scattered field Ug can also be separated by the 
Fourier transform in x, 



Plugging this representation back into ()2.1|) . and solving for each frequency ^ sep- 
arately, we find that u{ty) = a{^)e'^^^^'>y + 6(^)e-*''(«)''. The second term on the 
right corresponds to an incoming wave, which we insist must be zero since we want 
the scattered field to consist only of outgoing waves. Then Us{^,y) = a{£_)e^'^'^^'^y . It 
follows that Us{x,0) — J a(^)e'^^ d^, and 

z/?(eK(C,o)e'«"dC = (ru,)(x). 



The linear operator T (Dirichlet-to-Neumann map) then defines the relationship be- 
tween the traces Us|{y^o} a-iid (?j,Us|{y=o}- T{us\{y^o}) — {9y'Us)\{y=o}- On the bound- 
ary {y — 0}, the solution u = Ui + Ug should then satisfy 

dyU — Tu = dyUi — Tui + dyUg — Tus — g — Tf 
= 2g. 

Define the periodic domain (circle) 

S = R/27rZ. 

Define the first Brillouin zone K — [—^, To reduce the problem (|2.ip over to a 
family of problems over 5 x M, we define for g E L^(K^) the Floquet transform !F by 

Tig) = e-^"- " 2/)e''""", « e K- 

The sum can be considered as a Fourier series in the quasi- momentum variable a, 
with values in L^{S x R). The map g Tg is an isomorphism from L-^(R^) to the 
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direct product space L^{S x R) 9 . Floquet theory assures that the solution u can 
be written 

.1/2 

u{x,y) = / ?/„(a;,y)e*""da, (2.2) 

J -1/2 

where each function Ua is 27r-periodic in the x variable, and satisfies the equation 

AaUa + LO^ pUa = 0, 

where Aq = A + 2iadi — |a|^. The boundary condition dyU — Tu = 2g on {y = 0} 
translates to 

dyUa - TaUa = 2ga, On {y = 0}, 

where 

T^ua = i(3{n + a)wa(n)e^"'^, (2.3) 

(here Ua{n) are the Fourier series coefficients of Ua with respect to the s-variable) . 
Similar considerations apply at the lower boundary of the slab {y — —6}, where 
(assuming there is no incoming wave coming from below) we find dyUa + TaUa — 0. 

2.2. Existence and uniqueness of solutions. Let 51 = 5 x (0, — &), Tq = {y = 

0}, Th — {y = —6}. Define an admissible class of dielectric coefficients 

A = {p = pr+ ipt e L°°{rt) : Pro < Pr{x) < Pri and Pig < pi{x) < pi^ a.e.}, 

where pr^, Pig > 0. Given the incident wave Ui generated by the point source at (0, h), 
we must solve the family of problems 

Aq,Mq + uj'^PrUa + iuj^piUa = 0, in ri (2-4) 

(^ - Ta)Ua = 2ga, On To 

oy 

d 

(— +rQ)ua = 0, on Tf,, 

oy 

for all a G [— i, i]. Existence and uniqueness of weak solutions, with a uniform bound, 
may be obtained for pi > Q. 

Lemma 2.1. For each p ^ A with pi > and a G [—5,5],. problem {2.4\l admits a 
unique weak solution Ua € H^(fl). Furthermore, there exists a constant C depending 
on Pi, A, such that \\ua\\H'^(Q,) ^ independent of p £ A and a. 

Proof. For convenience we drop the subscript a on solutions. Define for u, v € 

a[u,v) — I Vu-\7v — uj'^ / puv + a^ / uv — 2ia / dxuv — I {Tau)v — I {Tau)v, 
Jn Jn Jn Jn Jfo 



and 



b{v) = 2 / gaV. 
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It is straightforward to show that a{u,v) defines a bounded sesquiHnear form over 
H^{n)xH^{^l), and that b{v) is a bounded linear functional on H^{n). Weak solutions 
u G H^{il) of (12. 4p solve the variational problem 

a{u,v)^b{v) for idlv eH^n). (2.5) 

The sesquilinear form a uniquely defines a linear operator A : H^(fl) — s- H^(fl) such 
that a{u,v) = {Au,v)H^{n), and the functional b{v) is uniquely identified with an 
element b € H^{il.) such that b{v) = {b,v) by reflexivity and an abuse of notation. 
Problem (|2.5p is then equivalently stated 

Au = b. (2.6) 

We intend to show that a is coercive by establishing a bound \a{u, u)\ > c > for 
aU u e H'^{Q) with ||u||_f/i(n) = 1- 

Integrating by parts in x, we see by periodicity that 

{dxu)u = 

Let A+{a) = {n e Z: /to(/3„) = 0} and A-{a) = Z - A+(a). Each A+(a) is a 
finite set and £ A+(a). We can write 




for j = 0, b. Notice that all /3„(a) S A (a) satisfy —i(3n > 0, so the second term on 
the right hand side of the equation above is real and non-negative. We then have 




Assuming — |Vup + |up — 1, and noticing that the first four terms 

on the right-hand side are purely real and the last two terms are purely imaginary, 
we find 




For convenience, write t = /^^ J2neA-(a) -i/3n(a)w(n)e™^M, r = /j-^(l -I- w^p^Dlw^, 
and s — Obviously t,r, and s are nonnegative real numbers which depend 

on u (and p in the case of r). Although t and s are essentially independent, r must 
satisfy 



(1 -I- pro)s < r < (1 -I- p,.Js. 



(2.7) 
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With this notation, 

2\a{u,u)\ > \1 + t + a^s — r| + uo^ piS. 

Note that in the case s > 2{i+p — J' have \a{u, u)\ > ^ui'^piS > ^^jfp ) ■ Otherwise, 
■s < 2(i+p — ) so that < |, and |a(u,M)| > i|l + f — r| > i. Hence, for all > 0, 
and all r satisfying (|2.7p . 

|aKu)|>c = min|^^^^,i|. 

The bound thus holds for every u with ||m||//i(o) = 1 and for every p £ A with 
Pi > 0. Given this coercivity bound, direct application of the Lax-Milgram Theorem 
yields existence of the bounded solution operator A^^ for problem (|2.6p such that 
< 1/c. Thus ||u||hi(o) < \\b\\HHn)/c. 

Given the bound on ||w||iji(f2), a uniform H^{fl) bound follows easily, since AqU = 
—Lo^pu is uniformly bounded in L^(f2). □ 

The complete solution u to the original problem p.ll) can then be reconstructed 
from (221). 

3. Optimal design. The goal of the optimization is to make a perfect image of 
the incident field Ui — H^^^Lur) on the opposite side of the slab {y < —b}. A "mirror 
image" to the incident field converging at the point (0, — (6 + hi)) would look like 

h'^\uj\/x'^ + (y + 6 + hi)"^), where H^'' — Hq^^ is the conjugate Hankel function. 
Thus we want the trace 

uix, -b) ^ H^^\uo^x-^ + hl) = q{x). 
The Bloch representations of u and / allows us to see that by setting 

Ua{x,-b) ^ qa{x), (3.1) 

with Qa defined similarly to fa, we get u{x, —b) = q{x). 

Problem (|2.4p together with the additional boundary condition (|3.ip is overposed. 
However, by allowing p to vary as a design variable, it may be possible to make (j3.ip 
hold approximately for each a. 

3.1. Problem Formulation. Let F{p,a) — Ualr^, where Ua £ H^{il) is the 
weak solution to problem (j2.4|) . 
Consider the minimization 

infj(p) = -/ \\Fip,a)-qJlda. (3.2) 

Theorem 3.1. The optimization problem has a solution. 

Proof. The proof follows the well-known direct method in the calculus of vara- 
tions. A is weak -k L°° compact. Consider a minimizing sequence {pn} with some 
subsequence (still denoted by {pn}) converging weak ★ to some p e A. De- 

note by Un the solution to the boundary value problem corresponding to pn- By 
Lemma (12.11) . the sequence {u„} has has a weakly convergent subsequence (still de- 
noted {un}), Un ^ u in H^{fl) for some u G H^{rt). Hence, u„ — > u strongly 
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p = l 




p = l 



unknown media 




h 

r=o 



r=-b 



Fig. 3.1. Model problem. A time-harmonic wave from a point source is incident from above. 
We wish to determine the unknown periodic medium such that a focus is attained below. 



in We hold v € H^{ft) fixed, and we have UnV — > uv strongly in L^{n). 



For k fixed, ap^{un,v) 



Qpj, {u, v) as n 



oo. Here we used the fact that 



rpa . iji/2(rj.) ^ if-V2(r^.) is continuous. Since uv e ap^{u,v) — > ap{u,v) 

as pk — > p weak ★ L°°. The trace map extends uniquely to a continuous linear 
operator r : H^{'^) — > if^/^(r). Thus, the traces are also convergent: Unlr^ ^ u\rj 
weakly in iJ^/^(r). This implies F{pn,a) F{p,a) weakly in iJ^/^(r), and, hence, 
Fa : A — > H^^^{r) is weak ★ L°° continuous for fixed a. But our bound on the 
solution is independent of a, and thus, the minimization problem has at least one 
solution p G A. Q 

3.2. Adjoint-state derivatives. Let 6p = Spr + iSpi be a "small" perturbation 
to the coefficient p. We denote the linearization of J(p) with respect to 5p by DJ{p)5p. 
We have 



DJ{p) 



For fixed a, 



Re{DFaip)Sp, F{p, a) - qa)L'^{Ti) da 

2 

^ Re{6p,DF*{p){F{a,a) - qc))L^n) da. 

2 



and 



DF*{p) : L^n) 
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DFa{p){5p) = dua where Su^ solves the linearized problem: 

AaSua + pSua — —Spu, in fl (3-3) 

d 

{- Ta)5ua = 0, on To 

dy 

d 

(— + Tq)5uq = 0, on Tf,, 

dy 

The adjoint of the derivative DFa{p){-) is the linear operator DF*{-) such 

that 

{DF^{p){5p),^)l2(t,) = {5p,DF*{p){i;))mn)- 
For all ip e L'^i^b) let Wa solve 



AaWa + OJ^ PrWa ~ ioj'^PiWa — 0, in H. (3-4) 

ii- - T*)wc. = 0, on To 

dy 

dy 

where T*f — —J2 */3„/(f^)e™'^. An integration by parts argument shows that 



SUa — / SpUaWa- 

Tt Jn 

The "gradient" of the functional J{p) is the function G{p) G L^(r2) for which 
DJ{p){5p) — Re ( / 6pG{p)] and G{p) — / UaWada 



where Wa solves p.4p with ^(J — F{p, a) ~ q^. 

Define the normal cones at the minimizer p = p,, + 'i-Pi- 



N{pr) = I? e [L^y : JJ^Pr - Pr)d^ > ypr^A), 

and 

N{p{) - |c e : jj^P^ - P^)dC > Vp, e ^ 

The gradient is normal to A at p: — Gp|-/V(p) 7^ 0. Thus, we have 

- / (Pr - Pr)Re{G{p))dx > Vp,. e A. 
Jn 

The continuity of the gradient is sufficient to reduce the above to pointwise optimality 
conditions. In particular, 

Pr = Pro =^ Re{G{p)) > 
Pro < Pr < Pn Re{G{p)) = 

Pr = Pn Re{G{p)) < 
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for almost every x in O. 
Similarly, 

/ (p, ~ p^)Im{Gip))dx > Vp, e A, 
Jn 

and the pointwise optimality conditions are: 

Pi = P»o => Im,{G{p)) < 
Pio < Pi < Pii Im{G{p)) ^ 

p, ^ p^, Im{G{p)) > 

for almost every x in 57. 

In a similar optimal design problem involving waveguides, one can use the weak 
continuation property to show [3J that at an optimal p, at least one of the four equality 
bound constraints above on the real and imaginary parts of p is attained at almost 
every point x E Q. A similar argument for the present problem is greatly complicated 
by the form of V J, involving an integral of a family of PDE solutions, rather than 
just a single background-adjoint pair. A precise a-priori characterization of optimal 
solutions is thus difficult. 

4. Numerical Results. Approximate solutions of (|3.2p are sought through nu- 
merical discretization and optimization. The variational problem (|2.5p was discretized 
with a first-order finite element method, using piecewise bilinear elements on a uni- 
form, rectangular grid. The design variable p was approximated by a piecewise con- 
stant function on the same uniform grid. The nonlocal boundary operators defined 
by (|2.3p were approximated by explicitly calculating the Fourier coefficients of the 
traces of the finite element basis, then truncating the sum in (|2.3|1 . The resulting 
finite element scheme can be shown to converge and to conserve energy, provided 
all the propagating terms are included in the sum This discretization leads to a 
large, sparse (except for the boundary terms), non-Hermitian matrix problem, which 
for simplicity is solved using the direct sparse solver in Matlab. 

The integral in (|3.2p was approximated by a discrete sum in a. By imposing 
X-axis symmetry in the designs, it suffices to integrate only over positive a. In the 
following examples, we used 20 equally-spaced positive values of a to approximate 
the integral. 

Despite the convenience of imposing a positive lower bound on the imaginary part 
of p in Lemma 2.1 for obtaining a uniform upper bound on solutions, we found that 
the numerical experiments were quite insensitive to small dissipations. Thus in most 
of the examples below, we set pig — 0. 

After discretizing J{p) through finite elements, optimization was accomplished 
with a straightforward projected gradient descent algorithm as in ,3,, using the adjoint 
as derived in Section 3 to calculate the gradient. We performed a large number of 
numerical experiments with this method, using different initial guesses for the design 
variable p, and varying the frequency lu, source and focus locations h and /ii, and 
constraints prg , Pn , Pig i Pii on the real and imaginary parts of p. 

Generally speaking, we found that the method was able to produce, from almost 
any initial guess, a structure which produced relatively high field intensity near the 
desired focus. Some parameter choices and initial guesses resulted in structures with 
much better focusing properties than others. 
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In the first experiment we start with a purely real material and allow the real 
part of the dielectric coefficient to vary between p^o — 1 ^^'^ Pri =12. The source is 
positioned at 2.5 units from the slab {h = 2.5), and we are looking to obtain a focus 
2.5 units on the opposite side of the slab (hi = 2.5). Through numerical optimization 
we discover the structure shown in Figure [47l1 which gives a spot size 0.424A. True 
subwavelength imaging is only possible if evanescent modes are present at the interface 
between the structure and the transmission medium. Figure [3] shows the evanescent 
modes for the generated image versus those of the target, clearly showing that the 
structure produces evanescent modes which approximate those of the objective. 

In the second example we image a source far away from the lens {h = 90), and we 
are looking for the lens that will produce the best image ten units from the lens {hi = 
10). Far away objects arc much more difficult to image, but through optimization of 
the structure, we are able to obtain a spot of size 0.38A. In the third example we start 
with a photonic crystal and allow purely real structures varying between p^o — 1 and 
Pri = 12. The source is positioned at 2.4 units from the slab {h = 2.4) and we are 
looking to obtain a focus four units away on the opposite side of the slab {hi = 4). 
The optimized structure is shown in Figure 14. 4[ and it gives a focus with spot size 
0.395A, which is much better than the one obtained if we just use photonic crystal as 
suggested by Luo et al. which produced a focus with spot size 0.67A |llj . 

Example four allows both the real and imaginary parts of the dielectric coefficient 
to vary (/o^o = li Pri = 12, and pi^ = 0, pi-^ = 1). The distance between the source 
and the lens and the distance between the lens and the image are set to four units 
{h = 4, hi ~ 4). The optimized structure gives a focus with a spot size 0.284A (Figure 
14. 5p . which is a significant improvement to those obtained by structures described in 
the research literature so far, although as far as we know, no real materials exist with 
these dielectric coefficients. 

All of the numerical experiments were computationally intensive. Each iteration 
required the solution of a family of diffi'action problems, two for each a, and because 
of the crude optimization method employed, many iterations were typically required. 
Most of the examples below took on the order of two days to run on a workstation. 
Since our purpose here was simply to illustrate the feasibility of designing such struc- 
tures through mathematical optimization, we did not devote much effort to improving 
the efficiency of the numerical methods. We believe that computation time could be 
improved by at least an order of magnitude with existing, but more sophisticated, 
numerical methods. 

5. Conclusions. We have demonstrated the feasibility of designing periodic 
structures with subwavelength focusing properties, via mathematical optimization. 
The approach is mathematically sound and through numerical discretization yields 
plausible, if somewhat non-intuitive, solutions. 

We point out two weaknesses with the approach presented here, all of which 
we believe could be improved through further work. First, the performance of the 
optimized structures tends to be very sensitive to small perturbations in material 
parameters. This fact combined with the relative complexity of the solutions means 
that attempting to fabricate these structures is at this point not an attractive idea. 
Both the sensitivity and the complexity could be addressed through adding appro- 
priate constraints or penalties to the objective. For example, optimization through a 
level-set approach would necessarily yield structures composed of only two materials, 
with no intermediate- index areas [TH]. This together with total variation penalties 
may yield "simpler" solutions. Of course such constraints may incur a decrease in the 
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Fig. 4.1. Results from Example 1. Upper left: initial (real) p; upper right: optimized solution; 
lower left: intensity cross section (blue) versus point source (red). Spot size is 0.424; lower right: 
real part of E field within the solution box. Eight periods are shown. 
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Fig. 4.2. Magnitude of evanescent modes for Example 1 (blue), versus those of the target (red). 
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Fig. 4.3. Results from Example 2. Upper left: initial (real) p; upper right: optimized solution; 
lower left: intensity cross section. Spot size is 0.380; lower right: real part of E field within the 
solution box. Eight periods are shown. 



performance of the structure. 

Second, the designs created with this approach are not translationaUy invariant. 
SpecificaUy, moving the structure laterally relative to the point source may result in 
decrease of focus. There are a few ways to circumvent this problem. In the simplest 
case, the distance h from the point source to the structure is large so that the wavefront 
impinging on the structure is nearly planar. Our experiments have shown that in this 
case focusing is not dependent on the lateral position of the point source (and is much 
less senstive to vertical translations), although the focus does translate periodically 
with the structure. One can also optimize for structures whose period is very small 
compared to the wavelength, although such structures tend to have less focusing 
power. The ultimate solution to this problem would require building the translation 
invariance of the focus (not of the medium) into the objective function. Unfortunately 
this would increase the complexity of the computations beyond the capability of the 
simplistic approach presented here, but with further work and refinement it should 
be possible. 
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lower left: intensity cross section (blue) versus point source (red). Spot size is 0.284; lower right: 
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